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Abstract 

We show how to apply a generalization of Local control design to 
the problem of saturation of a spin 1/2 particle by magnetic fields in 
Nuclear Magnetic Resonance. The generalization of local or Lyapunov 
control arises from the fact that the derivative of the Lyapunov function 
does not depend explicitly on the control field. The second derivative is 
used to determine the local control field. We compare the efficiency of 
this approach with respect to the time-optimal solution which has been 
recently derived using geometric methods. 



1 Introduction 

The control of spin systems by magnetic fields has a great potential for appli- 
cations in Nuclear Magnetic Resonance (NMR) [lj. An example is given by the 
preparation of a well-defined initial state such as the saturation process which 
consists in vanishing the magnetization vector of a sample by using an adequate 
pulse sequence [5J [3]. In the setting of NMR spectroscopy and imaging, the 
saturation control can be used for solvent suppression or contrast enhancement. 
Different approaches have been proposed up to date to solve this control prob- 
lem. They extend from intuitive schemes such as the Inversion Recovery (IR) 
method [1 to more elaborate pulse sequences based on geometric optimal con- 
trol theory [lj or optimization schemes [5J O HH [T7] . In a recent letter [3] , we 
showed that a gain of 60 % in the duration of the saturation process can be 
reached when using the time-optimal solution with respect to the IR one. In 
general, however, analytic solutions as derived in [3] are not available, and easily 
accessible numerical methods are required for the control of complex and realis- 
tic systems. Here, we investigate such a method and demonstrate its efficiency 
through a comparison with time-optimal analytic solutions [4J . Roughly speak- 
ing, the present method is an extension of local [7] or Lyapunov control [8], in 
terms of a so-called Lyapunov function V over the state space which is minimum 
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(or maximum) for the target state. In original Lyapunov control theory, a con- 
trol field that ensures the monotonic decrease of V is constructed with the help 
of the first time derivative V of V. Under well established general conditions on 
the Hamiltonian (see e.g. [5] for a recent review), the system converges asymp- 
totically towards the target state. This approach has been largely explored in 
quantum mechanics both from the physical 110] and mathematical points 
of view [T2] . 

In this paper, we will employ a target functional whose first time derivative 
does not depend on the available control fields. In this case, an inspection of V 
does not help to identify suitable control fields and standard local control theory 
cannot be used. Therefore, one has to resort to the second time derivative V, 
which allows the design of control pulses that speed up the increase of V [13] . 
On the first sight, one might assume that loosing control of V should be a 
disadvantage. However, it turns out that it also yields significant advantages: 
since control fields typically induce unitary dynamics, the independence of V 
on the control field implies the invariance of V under the corresponding unitary 
group. Consequently, the function V is constant for a continuous manifold of 
quantum states within which the control can induce dynamics. By not enforcing 
a given path through state space, but leaving the systems the opportunity to find 
its optimal path with a direction restricted only by the manifolds of constant 
V permits to not exert unnecessary control, but to leave the system substantial 
freedom to find its optimal path. Also, V allows to predict the dynamics of V 
for longer intervals than V, since it extrapolates for two infinitesimal time steps 
rather than one. This reduces the time-local character of Lyapunov control, 
and the efficiency of this generalized approach will be illustrated by different 
numerical simulations in this work. 

The paper is organized as follows. Section [2] deals with the description of 
the system and of the control problem. The description of the generalized local 
control approach and of the time-optimal one is done in Sec. [3J Different 
numerical results are finally presented in Sec. 2] for relevant situations in NMR. 
Conclusion and prospective views are given in Sec. [Sj 



2 The model system 

The state of a spin 1/2 particle in an NMR experiment is most conveniently 
expressed by its Bloch vector with components M x , M y and M z . If the radio- 
frequency magnetic field is in resonance with the spin, the Bloch equations are 
given by: 

dt ~ u V IVI z T 2 

^ = -u x M z -^ (1) 
= -cj y M x +uj x M y -M^Mo 

where oj x and uj v are the two components of the control field uJ, T\ and Ti are 
two relaxation rates characterizing the incoherent process due to the interaction 
with the environment and Mq is the thermal equilibrium magnetization. Since 
the maximum amplitude LJ max that the control fields uj can adopt induces a 
natural time-scale, it is convenient to introduce the scaled variables 

u = 2nuj/uj max , 
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r = 


{u m ax/2ir)t, 


r = 


2n/ {uJ max T 2 


7 = 


2n/ (cOmaxTl 


x = 


M/M . 



Since, in addition, the equations of motions are invariant under rotations around 
the z-axis, we can restrict the entire dynamics to the (y, z)-plane by assuming 
that uj y = 0. Doing so, Eqs. ([1} reduce to 

V = -Ty-uz 
z = 7(1 — z) + uy 

in terms of the scaled and scalar control field u = 2Ttu) x /ui max . 



3 Optimization schemes 

In the following we will take the thermal equilibrium state with M z = Mq and 
M x = My — as initial state. In terms of the scaled magnetization, this implies 
that the initial state is located at the north pole of the scaled Bloch sphere. As 
mentioned below, the goal of the control is to bring the scaled Bloch vector to 
the center of the sphere in minimum time with the constraint \u\ < 2ir on the 
control field. 



3.1 Generalized local control theory 

Since our goal is to reach a vanishing magnetization vector, which is equivalent 
to a maximum entropy state, we can choose any entropy-like functional as target. 
We consider the linear entropy defined as 

^(g)-l-Tr e 2 = i(l-x 2 -t/ 2 -z 2 ). (3) 

Reaching a maximum entropy in minimum time, i.e. as quickly as possible, can 
be done by choosing u such that the temporal increment of Si is maximized. As 
anticipated above, the time derivative 

^ = -iy + 7 (i- 2 )z, (4) 

of Si is independent of it, which is a direct consequence of the invariance of Si 
under coherent dynamics, as generated by u. We, therefore, have to inspect the 
curvature of Si , which essentially contains the information of how the impact of 
relaxation changes under the coherent dynamics. In the present case, it reads 



d 2 Si 
dt 2 



= uy{ 1 -2{ 1 -Y)z)+S^ (5) 



where S° denotes the curvature in the absence of a control field. Since maxi- 
mization of the curvature results in maximally fast increase of Si , which in turn 
implies an acceleration of the increase of entropy, we will choose, at any instance 
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of time, u such that it maximizes the curvature. Given that the modulus of u is 
bounded by 27r, this leads to the choice u — 2tt if \i > 0, and u = —2tt if \i < 0. 
As long as \l is finite, this implies that u is constant over a finite period of time 
so that the equations of motion can be integrated straightforwardly. Care is, 
however, necessary if /i = 0. In this case, one has to distinguish between stable 
solutions, where the control fields ensure that the condition fi = is preserved 
and unstable conditions, where the control, or the natural dynamics ensures 
that the condition // = is satisfied only for an isolated point in time. The 
latter case is realized for y = and it corresponds to u — {i.e. no control 
pulse), the system remaining in the line of equation y — to follow its natural 
dynamics. The former case is realized for z = z$ = 7/(2(7 — T)) and requires 
some more care. Let us, therefore, step back for a moment from the framework 
of differential equations and discuss Eqs. @ in terms of finite time steps. If we 
start out with fi < than the value of /1 will increase in time since we were as- 
suming the condition /1 = to be stable. Due to the finite time steps, however, 
we might end up in a situation where fi > 0. Now, in turn, /j, will decrease in 
time, so that we risk to end up in rapid changes of sign of //, which results in 
a rapidly changing sign of the control pulses. Decreasing the size of the time 
steps will decrease the size of the fluctuations of fj, and increase the frequency 
with which the control pulse changes sign. Once this frequency has exceeded all 
system frequencies sufficiently, we can replace the rapidly oscillating pulse by a 
time-averaged pulse 

u{t) = / di u(t)f(t - t) , (6) 

where f(t) is some function, like a Gaussian, that is positive in some finite 
domain around t = 0, and that decays sufficiently fast outside this domain. In 
the limit of infinitesimal time steps, the fluctuation of [i will vanish, and u is 
the pulse that has to be applied to preserve the condition /j, = 0. The control u 
needs to be chosen such that 



z = 7(1 - z) +uy = . (7) 
This implies that the local pulse solution can be written as 

u = - 1 (l-z) = -^(l-—t—) . (8) 

y y 2(7 - r) 

Note that u — > ±00 when y — > which defines a limit of admissibility due to 
the bound on the control given by 

\y\ > 17(7 - 2r)l (9) 
m - 27r(2r-2 7 )- iyj 

This means that for smaller values of y, the system cannot follow this horizontal 
line and a constant control u = ±27r has to be used. Together with Eq. ([2]) , the 
equation of motion for the y-component with u = u reads consequently as 

t = -Ty + - IS ' (10) 

y 4 (7 -r) 
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which has the solution 



ys = ± ] /e-^ -^) + ^ (11) 

where j/o is the initial y- coordinate on the horizontal axis. One deduces that u 

is determined through Eq. ([5]). and, u in turn, defines z(t) uniquely. One finally 
arrives at: 

u = 2ir for ^ > 0, (12) 

u = -2vr for fi < 0, (13) 

u = for y = and ^ = 0, (14) 

u = u as given by Eqs. ([5]) for /i = and z = zq, (15) 

and the complete solution is a concatenation of these pieces, as sketched in 
Fig. [TJ In the following, the control fields of maximum intensity will be called 
regular or bang, while those such that \u\ < 2tt will be said to be singular. 
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Figure 1: Temporal increment [y, i] of the Bloch vector in the y-z-plane in the 
presence of local control as given in Eqs. (|12|) - (fT5|) . The dashed arrow depicts 
the evolution along the line /i = according to Eq. (fT5]k the dotted arrow 
corresponds to free evolution according to Eq. (fl"4")h The curved, solid arrows 
correspond to Eqs. (fT2l and (flU)) . 

Using this material, we can follow the dynamics starting from the thermal 
equilibrium (y = 0, z = 1). This point is an unstable fixed point of the dynamics, 
i.e. following the above prescription one should not add any control pulse, the 
state being stationary. Any small deviation, however, will indicate that a finite 
control field needs to be applied. We will, therefore, begin with an infinitesimally 
displaced initial condition y(0) = — e, z(0) = 1. Then, according to Eq. (fT3"|) 
the system is driven with maximum intensity u = — 2n until the Bloch vector 
reaches a point with fj, = 0. At this point, the control field has to be changed 
such as to satisfy Eq. (fT5|) and the Bloch vector moves along the line fi = 
until its y-component vanishes. Once, this is achieved, the spin will relax to the 
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target state without any control field according to Eq. (fH|) . For the moment, 
we have been assuming that the control field can be chosen strong enough to 
ensure that the spin remains on the horizontal ji = 0- trajectory. Here, it is not 
the case due to the limitation in the field amplitude, a constant pulse according 
to Eq. (jT3J) has therefore to be applied when u — ±2ir along the horizontal line. 



3.2 Optimal control theory 

In this section, we briefly recall the tools of optimal control theory based on 
the application of the Pontryagin Maximum Principle (PMP). The reader is 
referred to |16| for a general application to dissipative quantum systems and 
to [H [13 \T7\ for examples in NMR. The first step of this method consists in 
introducing the pseudo-Hamiltonian % defined by 

H = P ■ (F (X) + uFtiX)) (16) 

where X = (y, z) and the vectors fields Fq and F\ have respectively the com- 
ponents (— r?/,7(l — z)) and (—z,y). The vector P is the adjoint state of coor- 
dinates (p y ,p z ). The PMP states that the optimal trajectories are solutions of 
the system 

X = Jl(X,P,v);P=-J±(X,P, V ) 
H{X,P,v) = max H(X,P,u) 

\u\<2tt 

H{X,P,v) > 0. 

Introducing the switching function $ = P ■ F\ , the PMP leads to two types of 
situations, the regular and the singular ones. In the regular case, $(X,P) ^ 
or vanishes in an isolated point while in the singular case $ is zero on a time 
interval. In the former situation, the corresponding control field is given by 
u = 2tt x sign[$]. In the latter case where <f> = (j> = cf> = ■ ■ • = 0, we introduce 
the vector V = (—7 + 72 — Tz, (-T + 7)2/) which satisfies <j> = P ■ V. The 
singular trajectories belong to the set of points where V is parallel to F\, i.e. 
to the points such that det(Fi, V) — 0. A straightforward computation leads to 
the two lines of equation y = and z = zq [3], which shows that the singular 
set is the same in the local control and optimal approaches. The corresponding 
optimal singular control is by definition the field such that the dynamics remains 
on the singular lines. Still this field is identical to local singular control. At this 
point, it is important to note that the trajectories given by the maximization 
of % are only extremal solutions, i.e. candidates to be optimal. Other tools like 
those introduced in [18] have to be used to get optimality results and to select 
among extremal trajectories the ones which are effectively optimal. 



4 Comparison of local and optimal control se- 
quences 

As a specific example to illustrate the efficiency of these two optimization 
schemes, we consider the control problem of a spin 1/2 particle coupled to a 
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thermal environment as analyzed and experimentally implemented in Ref. [I]. 
In this section, we will consider two examples characterized by different relax- 
ation parameters. In the first case, we have T2 = 8.8 ms and T\ = 61.9 ms 
which leads, with a maximum amplitude uj max = 2tt x 32.3 Hz, to T = 3.5 and 
7 = 0.5. In the second situation, we choose T2 — 6.2 ms and T\ = 61.9 ms 
which gives with the same maximum amplitude r = 5 and 7 = 0.5. 

Let us compare the time-optimal solution as computed in [?J I16j with the so- 
lutions obtained with the generalized local control theory. The numerical results 
are represented in Fig. [5] and [3] for the first and second cases, respectively. In 
the two cases, the time optimal solution is composed of two bang pulses of max- 
imum amplitude 2tt and two singular controls: starting from the initial point of 
coordinates (y = 0,z = 1), the first bang pulse drives the spin up to the horizon- 
tal singular line where /1 vanishes. The consecutive singular control ensures that 
the z-component remains constant until y — y Cl when the second bang pulse is 
applied. The position y c is computed to ensure that the switching function cf> 
and its first derivative <fi remain continuous along the optimal sequence. The 
values of y c are equal to -0.1605 and -0.1356 in the first and second cases. The 
second bang pulse increases the z-component of the magnetization vector un- 
til its y-component vanishes. The final part of the optimal sequence is a zero 
control along the vertical singular axis y = which allows to reach the target 
state. With the chosen relaxation parameters 7 and T, the time-optimal se- 
quence has a duration of 0.95050 and 0.97 in dimensionless units [defined above 
in Eq. $2$] in the first and second cases. When there is no bound on the control, 
the second bang does not exist and the optimal sequence is formed by a bang arc 
followed by the two singular trajectories along the horizontal and vertical lines. 
The solution of the generalized local control is very similar. Exactly like the 
time-optimal solution, it begins with a bang pulse followed by a singular pulse 
ensuring that fi vanishes, but this pulse is maintained longer than in the time- 
optimal solution. Only at y = —0.0862 and y = —0.0840 where \u\ — 2ir, the 
limited control intensity does not permit to maintain the condition ijl = and 
a bang pulse has to be applied. In the unbounded situation, note that the local 
and optimal solutions are identical. In the bounded case, if the second bang arc 
intersects the vertical z- axis, the dynamics relaxes without any control until 
its Bloch vector vanishes. In the case of Fig. [3J this leads to an overall control 
duration of 0.95098 which is slightly larger than the one of the time-optimal so- 
lution, the difference is of the order of 4.81 x 10~ 4 . Whereas local control does 
not reproduce the time-optimal solution exactly, we can conclude in this case 
that the longer duration of the former is negligible for all practical purposes. 
This conclusion is not true for all the values of the relaxation parameters. In 
particular in Fig. [31 since there is no intersection between the second bang pulse 
and the z- axis, we observe that the local control approach does not manage to 
reach the target, the minimum distance being equal to 9.9 x 10 -4 . 

Figure [4] gives a third example where the initial point of the dynamics is the 
south pole of the Bloch sphere, the relaxation rates being the ones of the first 
example. If we assume that the initial point is exactly on the z- axis then the 
local solution consists of letting act the relaxation with u — up to the center 
of the Bloch ball. Since this axis is unstable, any small deviation from this point 
will induces a bang control of amplitude ±2n up to the horizontal singular line. 
A second bang pulse is applied when the limit of admissibility of the control 
is reached and a zero control finally used up to the target state. The optimal 
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control sequence is quite similar, except for the second bang which is applied 
earlier than in the local control strategy. We observe in Fig. 0] that the local 
and the optimal solutions are very close to each other with a control duration 
respectively of 0.8758 and 0.8753, while if u — this time is of 1.3861. 
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Figure 2: (Color Online) (top): Evolution of the magnetization vector in the 
{y, z)-plane for the local control trajectory (blue or black) and the optimal 
solution (red or gray) with decoherence-times T2 = 8.8 ms and T\ = 61.9. 
The insert is a zoom of the trajectories near the singular horizontal line. The 
cross indicates the position of the limit of admissibility of the control along the 
horizontal axis. The bottom panel displays the corresponding control fields. 



5 Conclusion and perspectives 

The close-to-optimal performance of the generalized time-local control pulses 
reflect the potential of this new approach which can efficiently replace the stan- 
dard local control design when this strategy cannot be used. Whereas we have 
seen that sufficiently strong control typically permits to reach the target, it can 
also be missed with weak control fields. This opens up the general question of 
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Figure 3: (Color Online) Evolution of the magnetization vector similarly to Fig. 
[5] but with decoherence-times T2 = 6.2 ms and T\ — 61.9. 
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Figure 4: (Color Online) Evolution of the magnetization vector with 
decoherence-times T2 = 6.2 ms and T\ = 61.9 as in Fig. ® but for an ini- 
tial state corresponding to the south pole of the Bloch sphere. 
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convergence in the present generalized framework of Lyapunov control. In stan- 
dard Lyapunov control theory there is a well-established theory of convergence, 
based, e.g., on the concept of the Lasalle invariant set on which the derivative 
of the Lyapunov function vanishes [T^]. In addition, necessary and sufficient 
conditions on the Hamiltonian of the system have been derived to ensure the 
converge towards the target state [5] • Given the success of these tools for stan- 
dard Lyapunov control theory it seems within reach to also develop rigorous 
concepts to predict convergence properties of the present extension. 
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